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Abstract 

We use generalized Gaussian quadratures for exponentials to develop a new ODE 
solver. Nodes and weights of these quadratures are computed for a given bandlimit 
c and user selected accuracy e, so that they integrate functions e tbx , for all \b\ < c, 
with accuracy e. Nodes of these quadratures do not concentrate excessively near 
the end points of an interval as those of the standard, polynomial-based Gaussian 
quadratures. Due to this property, the usual implicit Runge Kutta (IRK) collocation 
method may be used with a large number of nodes, as long as the method chosen 
for solving the nonlinear system of equations converges. We show that the resulting 
ODE solver is symplectic and demonstrate (numerically) that it is A-stable. We use 
this solver, dubbed Band-limited Collocation (BLC-IRK), in the problem of orbit 
determination. Since BLC-IRK minimizes the number of nodes needed to obtain 
the solution, in this problem we achieve speed close to that of explicit multistep 
methods. 



1 Introduction 



Current methods for solving ODEs, be that multistep or Runge-Kutta, are 
based on polynomial approximations of functions. On the other hand, both 
recent and classical results [2fT|23f21|ll2f22] indicate that in many situations 
band-limited functions provide a near optimal tool for numerical integration 
and interpolation of functions. Using band-limited approximations, we con- 
struct a new method for solving the initial value problem for the ordinary 

1 This research was partially supported by AFOSR grant FA9550-07-1-0135, NSF 
grant DMS-0612358, DOE/ORNL grants 4000038129 and DE-FG02-03ER25583. 



Preprint submitted to Elsevier 



August 17, 2012 



differential equations (ODEs) and demonstrate certain advantages of such ap- 
proach. As an example, we consider orbit determination problem yielding a 
practical application as well as a gauge to ascertain the performance of new 
algorithms. 

It is well-known that choosing between equally spaced and unequally spaced 
nodes on a specified time interval (i.e., choosing a multistep vs a collocation 
based Runge-Kutta method), results in significantly different properties of 
ODE solvers. When applied to solving a constant coefficient, linear system 
of ODEs (i.e., computing a matrix or scalar exponential, the so-called test 
problem), multistep schemes are A-stable only if their order does not exceed 2, 
the so-called Dahlquist barrier. In contrast, an A-stable implicit Runge-Kutta 
(IRK) scheme may be of arbitrary order. A class of A-stable IRK schemes 
uses the Gauss-Legendre quadrature nodes on each time interval; the order of 
such methods is 2u, where v is the number of nodes (see, e.g., [9]). Having 
{TZe(z) < 0, z G C} as the region of absolute stability, assures that growth 
and decay of numerical solutions exactly mimics that of the analytic solutions 
of the test problem. A proper use of an A-stable scheme implies that the choice 
of parameters involves only accuracy consideration. 

Another property of interest, that of preservation of volume in the phase space, 
identifies symplectic integrators. Symplectic integrators preserve a particular 
conserved quantity of nonlinear Hamiltonian systems as well as an approxi- 
mate Hamiltonian. In problems of orbit determination, a symplectic integrator 
would keep the correct orbit more or less indefinitely with the error accumu- 
lating only in a position along the accurate orbit, thus closely reproducing a 
particular behavior of analytic solutions of nonlinear Hamiltonian systems. We 
note that IRK schemes which use the Gauss-Legendre nodes are symplectic 
(see, e.g., |9]). 

While IRK schemes with the Gauss-Legendre nodes provide an excellent dis- 
cretization of a system of ODEs, using many such nodes on a specified time 
interval presents a difficulty. It is well known that the nodes of the Gauss- 
Legendre quadratures (as well as any other polynomial based Gaussian quadra- 
tures) accumulate rapidly towards the end points of an interval. A heuristic 
reason for such accumulation is that these quadratures have to account for a 
possible rapid growth of polynomials near the boundary. 

In this paper we demonstrate that, within IRK collocation methods, quadra- 
tures based on polynomials may be replaced by quadratures for band-limited 
exponentials. The nodes of these quadratures do not accumulate significantly 
toward the end points of an interval (a heuristic reason for an improved ar- 
rangement of nodes is that the exponentials do not grow anywhere in the 
interval). The quadratures for exponentials have been successfully used in 
problems of wave propagation [2] (see also [16]). We note that in problems of 
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wave propagation solutions are naturally well approximated by band-limited 
functions. It is also typical that solutions of ODEs are well approximated by 
band-limited exponentials. However, for ODEs there may be exceptions since 
solutions of some ODEs may, in fact, be polynomials or other functions that do 
not have an efficient approximation via exponentials and where polynomial- 
based quadratures are more effective. Thus, our method addresses numerical 
integration of ODEs whose solutions are well approximated by band-limited 
exponentials. 

Unlike the classical Gaussian quadratures for polynomials which integrate 
exactly a subspace of polynomials up to a fixed degree, the Gaussian type 
quadratures for exponentials use a finite set of nodes to integrate an infinite 
set of functions, namely, ie* b M on the interval \x\ < 1. As there is no 

l J |6|<c 

way to accomplish this exactly, these quadratures are constructed so that all 
exponentials with |6| < c are integrated with accuracy of at least e, where e 
is arbitrarily small but finite. Such quadratures were constructed in [1] and, 
via a different approach in [23J. As observed in j2], quadrature nodes of this 
type do not concentrate excessively toward the end of the interval; the density 
of nodes increases toward the end points of the interval only by a factor that 
depends on the desired accuracy e but not on the overall number of nodes. 

Using quadratures integrating with accuracy e 2 band-limited exponentials with 
banlimit 2c, we naturally arrive at a method for interpolation with accuracy 
e of functions with bandlimit c (see [2"3"1fT] ). It turns out that the nodes and 
weights of quadratures with interpolation accuracy e w 10~ 15 can be con- 
structed using only the standard double precision machine arithmetic. How- 
ever, computing the integration matrix for accuracy e ~ 10 -15 does require an 
extended precision. We describe an approach that uses only quadruple preci- 
sion for this purpose. Importantly, we use only the standard double precision 
when applying such quadratures and integration matrices within BLC-IRK 
method. 

While analytically the classical Gauss-Legendre quadratures for polynomials 
are exact, in practice their accuracy is limited by the machine precision. By 
choosing (interpolation) accuracy e ~ 10~ 15 , our integrator is effectively "ex- 
act" within the double precision of machine arithmetic. Remarkably, using 
a particular construction of the integration matrix, we show that BLC-IRK 
method is (exactly) symplectic and, with high accuracy, A-stable. This result 
was unexpected and indicates that properties of approximate quadratures for 
band-limited exponentials need to be explored further. 

While IRK schemes require solving a system of nonlinear equations at each 
time step, it does not automatically imply that such schemes are always com- 
putationally more expensive than explicit schemes. In the environment where 
the cost of function evaluation is high, the balance between the necessary it- 
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eration with fewer nodes of an implicit scheme vs significantly greater number 
of nodes of an explicit scheme (but no iteration), may tilt towards an implicit 
scheme. In problems of orbit determination, we use an additional observation 
that most iterations can be performed with an inexpensive (low fidelity) grav- 
ity model, making implicit schemes with a large number of nodes per time 
interval practical. We select the problem of orbit determination as an example 
where our approach is competitive with numerical schemes that are currently 
in use (see [3]). We take advantage of the reduced number of function calls 
to the full gravity model in a way that appears difficult to replicate using 
alternative schemes. 

In order to accelerate solving a system of nonlinear equations, we modify the 
scheme by explicitly exponentiating the linear part of the force term. For 
the problem of orbit determination, this modification effectively makes use of 
the fact that the system is of the second order to accelerate convergence of 
iterations. So far we did not study possible acceleration of iterations using 
spectral deferred correction approach as in |5fT4~f8][TU] . 

We start by providing background information on quadratures for band-limited 
functions in Section [2J We then describe BLC-IRK method in Section [3] (with 
some details deferred to Appendix). In Section H] we detail our algorithm and 
provide examples. 



2 Preliminaries: quadratures for band-limited functions 



2.1 Band-limited functions as a replacement of polynomials 



The quadratures constructed in [23,1] break with the conventional approach of 
using polynomials as the fundamental tool in analysis and computation. The 
approach based on polynomial approximations has a long tradition and leads 
to such notions as the order of convergence of numerical schemes, polynomial 
based interpolation, and so on. Recently, an alternative to polynomial approx- 
imations has been developed; it turns out that constructing quadratures for 
band-limited functions, e.g., exponentials e lbx , with \b\ < c, where c is the 
bandlimit, in many cases leads to significant improvement in performance of 
algorithms for interpolation, estimation and solving partial differential equa- 
tions PITMT] . 
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2.2 Bases for band-limited functions 



It is well-known that a function whose Fourier Transform has compact support 
can not have compact support itself (unless it is identically zero). On the other 
hand, in physics duration of all signals is finite and their frequency response 
for all practical purposes is also limited. Thus, it is important to identify 
classes of functions which are essentially time and frequency limited. Towards 
this end, it is natural to analyze an operator whose effect on a function is 
to truncate it both in the original and the Fourier domains. Indeed, this has 
been the topic of a series of seminal papers by Slepian, Landau and Pollak, 
|22|I12|I13|I18|I19|I20|I21| . where they observed (inter alia) that the eigenfunctions 
of such operator (see (j2J) below) are the Prolate Spheroidal Wave Functions 
(PSWFs) of classical Mathematical Physics. 

While periodic band-limited functions may be expanded into Fourier series, 
neither the Fourier series nor the Fourier integral may be used efficiently for 
non-periodic functions on intervals. This motivates us to consider a class of 
functions admitting a representation via functions of the form , x G 

[—1, 1], with a fixed parameter c (bandlimit). Following [lj, let us consider the 
linear space of functions 

£ c = jueL°°([-l,l]) | u(x) = J2a k e tcbkX :{a k } keZ el 1 , b k e[-l,l]\. 

Given a finite accuracy e, we represent the functions in £ c by a fixed set of 
exponentials {e tCTkX } k=1 , where M is as small as possible. It turns out that by 
finding quadrature nodes {r k } jf =1 and weights {w k } %L 1 for exponentials with 
bandlimit 2c and accuracy e 2 , we in fact obtain (with accuracy e) a basis for 
S c with bandlimit c pQ. 

The generalized Gaussian quadratures for exponentials are constructed in [1] 
(see [23J for a different construction), which we summarize as 

Lemma 1 For c > and any e > 0, there exist nodes — 1 < t\ < T2 < ■ ■ ■ < 
T/y/ < 1 and corresponding weights w k > 0, such that for any x G [—1,1], 



i M 
J i e ictx dt-J2 w k e lCTkX 



k=l 



<e, (1) 



where the number of nodes, M , is (nearly) optimal. The nodes and weights 
maintain the natural symmetry, r k = —tm-u+i and w k = WM-k+i- 

Remark 2 The construction in £iy is more general and yields quadratures 
for band-limited exponentials integrated with a weight function. If the weight 
function is 1 as in LemmaU\ then the approach in fTj identifies the nodes of 
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the generalized Gaussian quadratures in (TJP as zeros of the discrete prolate 
spheroidal sequences (DPSSs) JW$ , corresponding to small eigenvalues. 

Let us now consider band-limited functions, 

B c = {fe L 2 (R) | f(u) = for \co\ > c}, 

and briefly summarize some of the results in [2"2"fl2|ll3|ll8f21] . Let us define the 
operator F c : L 2 [-1, 1] L 2 [-1, 1], 

F c (^)(w) = j\ icx ^(x)dx, (2) 

where c > is the bandlimit. We also consider the operator Q c = ^F*F C , 
^ / , w s 1 Z" 1 sin(c(w — a;)) , , . , 

QcW(y) = - / — ^ -V> * ^- (3) 

The eigenf unctions ^(j, ^J, ' ' ' °f Qc coincide with those of F c , and the 
eigenvalues fij of Q c are related to the eigenvalues Xj of F c as 

^ = ^-|Aj-| 2 , J = 0,1,2,.... (4) 

While all /Xj < 1, j = 0, 1, . . . , for large c the first approximately 2c/tt eigen- 
values fij are close to 1. They are followed by O(logc) eigenvalues which decay 
exponentially fast forming a transition region; the rest of the eigenvalues \ij 
are very close to zero. 

The key result in [22] states that there exists a strictly increasing sequence of 
real numbers 70 < 71 ... , such that tpj are eigenfunctions of the differential 
operator, 

L c r 3 = (-(1 - x 2 ) ^ + 2xj- + c 2 * 2 ) r 3 {x) = ■ ( 5 ) 

The eigenfunctions of L c have been known as the angular prolate spheroidal 
functions (PSWF) before the connection with ([2D was discovered in [22] (by 
demonstrating that L c and F c commute). We note that if c — > 0, then it 
follows from fl5]) that, in this limit, ip c , become the Legendre polynomials. 
In many respects, PSWFs are strikingly similar to orthogonal polynomials; 
they are orthonormal, constitute a Chebychev system, and admit a version of 
Gaussian quadratures [23J. 

Since the space £ c is dense in B c (and vice versa) [lj, we note that the 
quadratures in [23j may potentially be used for the purposes of this paper 
as well (the nodes of the quadratures in |23j and those used in this paper 
are close but are not identical). Importantly, given accuracy e, the functions 
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?/>o, ipti ^21 ' ' ' i V'm-i ma y be used as a basis for interpolation on the interval 
[—1,1] with 7~i, r 2 , • • • , tji as the interpolation nodes, provided that these are 
quadrature nodes constructed for the bandlimit 2c and accuracy e 2 . Given 
functions ip^, ijj^, ■ ■ ■ , V'm-d we can construct an analogue of Lagrange in- 
terpolating polynomials, R%(x) = Djio 1 a kj^Pj(x), x G [—1, 1], by solving 

M-l 

5 kl = Rl{n) = £ o^fa) (6) 
i=o 

for the coefficients a^j. The matrix ipj(ri) in ([S]) is well conditioned. 

A well-known problem associated with the numerical use of orthogonal poly- 
nomials is concentration of their roots near the ends of the interval. Let us 
consider the ratio 

r(M,e) = , (7) 

T [M/2] — T LM/2J-1 

where "|_M/2J" denotes the least integer part, and look at it as a function of 
M. Observing that the distance between nodes of Gaussian quadratures for 
exponentials changes monotonically from the middle of an interval toward its 
end points, and that the smallest distance occurs between the nodes closest 
to the end point, the ratio may be used as a measure of node accumula- 
tion. For example, the distance between the nodes near the end points of the 
standard Gaussian quadratures for polynomials decreases as 0(1/M 2 ), so that 
we have r(M, e) = 0(1/M), where M is the number of nodes. In Figured] we 
illustrate the behavior of r(M, e) for the nodes of quadratures for band-limited 
exponentials. This ratio approaches a constant that depends on the accuracy 
e but does not depend on the number of nodes. 

Another important property of quadratures for exponentials emerges if we 
compare the critical sampling rate of a smooth periodic function, to that of 
smooth non-periodic function defined on an interval. Considering bandlimit 
c as a function of the number of nodes, M, and the desired accuracy e, we 
observe that the oversampling factor, 

<*{M,e = — — >1, 
c(M, e) 

approaches 1 for large M. We recall that in the case of the Gaussian quadra- 
tures for polynomials, this oversampling factor approaches | rather than 1 
(see e.g. [7]). 
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Figure 1. The ratio r(M, e) in (J7J) as a function of the number of nodes M and 
interpolation accuracy e ~ 10~ 35 (top curve, dashed), e ~ 10 -8 ' 5 (middle curve, 
dotted) and e ~ 1CP 13 (middle curve, solid). The dots on the solid curve indicate 
the number of nodes of quadratures used in our numerical experiments. The bottom 
curve shows this ratio for the Gauss-Legendre nodes. 

2.3 Interpolating bases for band-limited functions 

A basis of interpolating band-limited functions for the bandlimit c and accu- 
racy e plays the same role in the derivation of a system of nonlinear equations 
for solving ODEs as the bases of Lagrange interpolating polynomials denned 
on the Gauss-Legendre nodes. While ([6]) relies on available solutions of the 
differential equation ([5]), interpolating basis functions may also be obtained 
by solving the integral equation (T5]) (see [1]). 

We start by first constructing a quadrature for the bandlimit 2c > and 
accuracy threshold e 2 > 0, yielding M nodes {r m }^ =1 and weights {w m }^ =1 . 
For the inner product of two functions f,g £ S c , we have 



Following [T], we discretize (j2J) using nodes {r m } m=1 and weights {w m } 
and obtain an algebraic eigenvalue problem, 




M 





1=1 
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The approximate PSWFs on [—1, 1] are then defined consistent with (jSJ) as 



j M 



= - y £w l e iaH t9 j (<n), (9) 

'h i=i 

where rjj are the eigenvalues and Vf^Tj) the eigenvectors in (jSJ). Following [T], 
we then define the interpolating basis for band-limited functions as 



M 

Rk(x) = J2r k ie tCT ' x , fc = l,...,M, (10) 

where 

r k i = Y, w ^j(n)—^j(n)wi. (n) 

j=i 

It is shown in [1] that the functions Rk(x) are interpolating, Rk{Ti) = 8j~i. 



3 BLC-IRK method 

3. 1 Discretization of Picard integral equation 

We consider the initial value problem for a system of ODEs, 

y' = f(t,y), y(0)=y o , 

or, equivalently, 

y(t)=y + /*f(s,y(s)) ds. (12) 
Jo 

It is sufficient to discretize (fT2l) on the interval [0, t] since, by shifting the time 
variable, the initial condition may always be set at t = 0. We require 

y'(tr J ) = f(tr J ,y(tT,)),...j = l,...,M, 

where {rj}jL 1 are Gaussian nodes for band-limited exponentials on [0, 1] (con- 
structed for an appropriate bandlimit c and accuracy e). We approximate 

M 

||f(tr,y(*r))-^f(tr i ,y(tr i ))i? i (r)|| < e, r £ [0, 1] (13) 

i=l 

where Rj(r) are interpolating basis functions associated with these quadra- 
tures and briefly described in Section [2.31 (see [Tf2] for details). Using (TTSj) . 
we replace f in (TI2"| and evaluate y(tr) at the quadrature nodes yielding a 
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nonlinear system, 

M Ti M 

y(foi) = yo+J2 f ( tr i' y( tT j)) / Rj(s)ds = yo+Y, SijfitTj, y(tTj)), i = 1, . . . M, 
3=1 Jo 3=1 

(14) 

where S^- = /J* Rj(s)ds is the integration matrix. Solving for {yitTj)}^, we 
have from (j!2p 

M 

y(0 = yo + E w i f (^>y(^))> ( 15 ) 

3=1 

where {wj}jL 1 are the quadrature weights. The result is an implicit Runge- 
Kutta method (IRK) where the usual Gauss-Legendre quadratures are re- 
placed by Gaussian quadratures for band-limited exponentials. 

The nodes, weights, and the entries of the integration matrix are typically 
organized in the Butcher tableau, 



T 



s 



Unlike in the standard IRK method based on Gauss-Legendre quadratures, we 
solve (HM on a time interval containing a large number of quadrature nodes, 
since these nodes do not concentrate excessively near the end points. This 
implies that the interval [0, t] may be selected to be large in comparison with 
the usual choices in RK methods. 



3.2 Exact Linear Part 



In many problems (including that of orbit determination), the right hand side 
of the ODE, f(£, y), may be split into a linear and nonlinear parts, 

f(t,y(t)) = Ly(t) + g(t,y(f)), 

so that the integral equation ( 1T21) may be written as 

y(t) = e tL y + f e^ L g(s, y(s)) ds. (16) 
Jo 

If the operator e tL can be computed efficiently, this formulation leads to savings 
when solving the integral equation iteratively. 

We discretize ( jT5j) by using ( 1T31) and obtain 
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y(tr t ) = e tT ' L y + ^2e t{n -^ L g(tr j ,y(tT J )) / * R^ds 

3=1 J ° 

M 

= e^ L yo + £ ^-e^-^gfe y(tr,)) (17) 

3=1 

where S{j = Rj(s)ds. We note that (|14|) is a special case of f TTTj) with L = 
and g = f . 

J.J Symplectic integrators 

Following [IT], let us introduce matrix Ai = {^Ar/}fcj=i for an M-stage IRK 
scheme, 

rrikj = w k Skj + WjSjk - w k wj, (18) 
where the weights w = {w k }^ =1 and the integration matrix S = {S k j}jfj =1 
define the Butcher's tableau for the method. 

It is shown in [T7] that 

Theorem 3 If matrix M = flin (COD, then an M-stage IRK scheme is sym- 
plectic. 

This condition, M. = 0, is satisfied for the Gauss-Legendre RK methods, 
see e.g. |lfTT] . We enforce this condition for BLC-IRK method by an 0(e 2 ) 
modification of the weights and of the integration matrix. For convenience, we 
consider the band-limited exponentials and integration matrix on the interval 
[—1, 1] (rather than [0, 1] usually used for ODEs). 

Proposition 4 Let {jj} 1 ^ =l and {wj}jL 1 be quadrature nodes and weights 
for the bandlimit 2c and accuracy e 2 . Consider interpolating basis functions 
on these quadrature nodes, Ri(r), Ri{Tj) = 5ij, i,j = 1,...,M, and define 
Fi(r) = JZ l Ri(s) ds. Then we have 



1 M 

/ F^Flir) dT-J2w k F 3 (T k )F'(T k 
k=i 



< ^ 



(19) 



or 



and 



or 



(^J Rj(s) ds^j Ri(r) dr — w,i J Rj(s) ds 



< t 



1 M 

j Ri(r) w k Ri{r k 
J ~ 1 k=i 



(20) 



Ri(s) ds — uii 



< t 
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PROOF. The relations in (TT9]) and (120]) is the property of the quadrature, 
since the bandlimit of the product Fj(r)F/(r) is less or equal to 2c and that 
of Ri{r) is less or equal to c. Due to the interpolating property of Ri(r), we 
have 



M M 



J2 w k F,(T k )Fl(r k ) = E(r Rj(s) ds) wuR^n) = Wi P Rj(s) ds (21) 

J 1 7 1 V*'— 1 / J —1 



k=l k=l 

and 



M 

^w k Ri(r k ) = Wi 

k=l 

Also, by definition, 

f^F^FUr) dr = £ (jTi^(«) da) ^(r) dr, 
and the result follows. 

Theorem 5 Let {rj}jL l be quadrature nodes of the quadrature for the ban- 
dlimit 2c and accuracy e 2 and Ri{t), Ri(Tj) = 5ij ; i, j = 1, . . . , M , the corre- 
sponding interpolating basis. Let us define weights for the quadrature as 

w k = J_ R k {r)dr (22) 

and the integration matrix as 

w k 

Then 

WkSki + mSik - w k wi = 0, (24) 
and the implicit scheme using these nodes and weights is symplectic. 



PROOF. Using the propositions above, we observe that the weights defined 
in (|2"2"|) are the same (up to accuracy e 2 ) as those of the quadrature. The result 
follows by setting F k {r) = JLiR k (r) dr, F' k {r) = R k {r) and integrating by 
parts to obtain 

w k S ki + wiS [k - w k wi = J ^ Fi(t)FI.(t) dr + J ^ F fc (r)F/(r) dr - w k w t 

= F,(l)F fc (l) - w k wi. 

By the definition of the weights, we have F k (l) = w k and, hence, Fi(l)F k (l) — 
w k wi = 0. 
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3.4 Construction of the integration matrix 



There are at least three approaches to compute the integration matrix. Two of 
them, presented in Appendix, rely on Theorem [5] and differ in the construction 
of interpolating basis functions. In what appears to be a simpler approach, the 
integration matrix may be obtained without computing interpolating basis 
functions explicitly and, instead, using a collocation condition together with 
the symplectic condition (1241) . 

Besides satisfying ( 124")) . we require that the method accurately solves the test 
problem 



y' = icr k y, y(-l) = 

on the interval [—1, 1], where r k are the nodes of the quadrature, k — 1, . . . , M. 
To find solution of this equation, yk{t) = e lCTkt — e~ lCTk , we require (see (TL41) ) 
that at the nodes, t = r, 



e iCT k Tj _ e 



-icr k M 



%CT k 



,JCT k T m 



m=l 



<e, 



(25) 



and compute the integration matrix as the solution of ( 124)) satisfying an ap- 
proximate collocation condition ( 125)) . 



Defining the symmetric part of the integration matrix 



T 



jin 



WjW m 



(26) 



we set 



so that it follows from (124")) that Aj m is anti-symmetric, 



(27) 



Aj m A m j 0. 



From (125]) . we obtain 



e icr k Tj _ e ~icr k M 



M 



lCT k 



E T , 



lCT k T m 



jrn^m^ 



lCT k Tm 



m=l 



m=l 



(28) 



We split the real and imaginary parts 



^ . Tj m € Ujk iVjk) 



lCT k 
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where 

M 

U jk = (Tj + 1) Sine {cT k (Tj + l)/2) COS [cT k {j j - l)/2) - T j™ COS (cT k T m ) 

771=1 

and, since T jm = Tj(M- m +i)due to symmetry of the weights, 

Vjk = {Tj + 1) sine (cr k (Tj + l)/2) sin (ct^tj - l)/2) . 
Splitting the real and imaginary parts in (128]) . we have 

M M 
U jk = A jmW m COS (cT k T m ) , V jk = ^ A jmW m sin (cT k T m ) . 

771=1 777 = 1 

Since matrices cos (cT k T m ) and sin (cr k r m ) are rank deficient, we choose to 
combine these equations as follows 

M 

Ujk + V jk = A jm W m (COS (cT k T m ) + SU1 (cT k T m )) . (29) 

777=1 

Solving this equation (using quadruple precision since this system is ill-conditioned), 
we find matrix A and discover that, while Sj m = Tj m + Aj m w m makes f l25|) 
into an equality, it is not exactly anti-symmetric. We then anti-symmetrize 
this matrix by averaging, Aj m = —A m j = (^Aj m — A m j^j /2. Combining so ob- 
tained anti-symmetric matrix A with (126|) via (j27|) . we obtain the integration 
matrix. We then verify that the inequality (1251) is satisfied. 

Remark 6 The fact that integration matrix satisfies [24\ ) and the inequality 
l[2l)\) indicates that, perhaps by a slight modification of nodes and weights of 
the quadrature, it might be possible to satisfy (24\ ) and I[2d*\) with e = 0. 

3.5 A-stability of the BLC-IRK method 

As shown in e.g. O Section 4.3], in order to ascertain stability of an IRK 
method, it is sufficient to consider a rational function 

r( z ) = 1 + zw\l - zS)- l l, (30) 

where S is the integration matrix, w is a vector of weights and 1 is a vector 
with all entries set to 1, and verify that \r(z)\ < 1 in the left half of the 
complex plane, IZe (z) < 0. This function is an approximation of the solution 
e z of the test problem 

y' = zy, y(o) = i 

computed via (]T4"]) and (]T5]) on the interval [0, 1]. If all poles of r(z) have a 
positive real part, then it is sufficient to verify this inequality only on the 
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Figure 2. Eigenvalues (computed using quadruple precision) of the integration matrix 
for BLC-IRK scheme with 64 nodes corresponding to the bandlimit 17n and, for 
comparison, eigenvalues of the integration matrix of the standard IRK scheme 64 
Gauss-Legendre nodes. 

imaginary axis, z = iy, y G K. In fact, it may be possible to show that r(z) is 
unimodular on imaginary axis, \r(iy)\ = 1, for y G KL Implicit Runge-Kutta 
methods based on Gauss-Legendre nodes are A-stable (see e.g [9j) and, indeed, 
for these methods r(z) is unimodular on imaginary axis. 

Given an M x M matrix S with Mi complex eigenvalues and M 2 real eigen- 
values, if function ( 1301) is unimodular on the imaginary axis then it has a 
particular form, 

r(z) = n -^-^r n -^i, (3i) 

k=l 2 ~ A k Z - \ k k'=l Z _ A k' 

where 2M\ + M2 = M. Currently, we do not have an analytic proof of A- 
stability of BLC-IRK method; instead we verify this property numerically. We 
compute eigenvalues of the integration matrix to obtain the poles of r(z) and 
check that all eigenvalues have a positive real part separated from zero. For 
example, integration matrix for BLC-IRK method with 64 nodes (bandlimit 
c = 17n) has all eigenvalues with real part larger than 0.7 • 10~ 3 . One way 
to check that r(z) has the form ( IHTj) is to compute r(— A^, 1 ) and r(— A^ 1 ) to 
observe if these are its zeros. In fact, it is the case with high (quadruple) 
precision. 



One can argue heuristically that since a rational function with M poles has 
at most 2M real parameters (since matrix 5* is real its eigenvalues appear in 
complex conjugate pairs) and since, by construction, r(iy) for \y\ < c is an 
accurate approximation to e iy (which is obviously unimodular), r(z) is then 
unimodular. It remains to show it rigorously; a possible proof may depend on 
demonstrating a conjecture in Remark El 
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4 Applications 

4-1 Algorithm 

A fixed point iteration to solve (TTTj) is the usual choice (as long as it converges). 
Let Nn denote the number of iterations, which can either be set to a fixed 
number or determined adaptively. Labeling the intermediate solutions in the 
iteration scheme as i = 1, . . . , N it , we have 

(1) Initialize y^fa) = y , i = 1, . . . , M. 

(2) For k = l,...,N it 
For m — 1, . . . , M 

(a) Update the solution at the node m: 



(b) Update the right hand side: g(£r m , y^ k '(tT m )) 

We note that the updated value of y^ k \tT m ) is used in the computation at 
the next node r m+i within the same iteration k. This is essential for a fast 
convergence. 

Remark 7 Although we currently apply the integration matrix directly, using 
a large time interval and, consequently, a large number of nodes per interval, 
opens a possibility of using fast algorithms for this purpose. Such algorithms 
may be faster than the direct application of the matrix only for a sufficiently 
large matrix size and are less efficient than the direct method if the size is 
relatively small. Since we may choose many nodes, it makes sense to ask if 
the integration matrix of an IRK type method may be applied in O (M) or 
O(MlogM) operations rather 0(M 2 ). We mention an example of an algo- 
rithm for this purpose using the partitioned low rank (PLR) representation 
(as it was described in e.g., f^) but leave open a possibility of more efficient 
approaches. 



4-2 Problem of Orbit Determination 

Let us consider the spherical harmonic model of a gravitational potential of 

degree N, 




e tr m L yo + e t(r m -^)L g(tT . ; y(*)(f Ti )) 




(32) 
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with 

n 

Y n {9,\) = P™{sm8){C nm cos{m\) + S nm sm{m\)), (33) 

m=0 

where P™ are normalized associated Legendre functions and C nm and S nm 
are normalized gravitational coefficients. In case of the Earth's gravitational 
model, n is the Earth's gravitational constant and R is chosen to be the Earth's 
equatorial radius. Choosing the Cartesian coordinates, we write ( r ),r = 
(x,y,z), assuming that the values V"W (r) are evaluated via f l3"2"j) by chang- 



ing from the Cartesian to the spherical coordinates, r = y/x 2 + y 2 + z 2 , 9 = 
arcsin(z/r) and A = arctan {y/x). 

We formulate the system of ODEs in the Cartesian coordinates and denote the 
solution as r(t) = (x(t),y(t), z(t)). Setting (r) = W (Af) (r), we consider 
the initial value problem 



dt 2 



r ( t ) = _ G W ( r (t)) , r (0) = r 



(34) 



We observe that in gravitational models the first few terms are typically large 
compared to the rest of the terms. For example, in EGM96 gravitational model 
|15| the only non-zero coefficients for Y 2 (9,X) are C20 C*22 and S22, where 
C 20 « -0.48 -10- 3 , C 22 w 0.24 -10" 5 , and ^ 22 ps -0.14- 10" 5 , whereas all other 
coefficients are less than 0.14 • 10~ 5 . For this reason it makes sense to split the 
force as 

G (A0 (r) = G (2) (r) + ^ G (N) (r) _ G (2) (r) ^ 

and use only (r) in most of the iterations (since using the full model, 
Q( N ) ( r ) ; may be expensive). 

We first use the gravity model of degree iV = 2 on a large portion of an orbit 
(e.g., 1/2 of a period) to solve the system of nonlinear equations via fixed 
point iteration. Once the approximate solution r(t) to 



d 2 ~, x 



-G( 2 )(f(t)), r(0) = r 



( x ^ 

Xq 

\ Zo J 



is obtained, we then access the full gravity model (r(tTj)) to evaluate the 
forces at the nodes Tj which, by now, are located close to their correct posi- 
tions. We continue iteration (without accessing the full gravity model again) 
to adjust the orbit. This results in an essentially correct trajectory. At this 
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point we may (and currently do) access the full gravity model one more 
time to evaluate the gravitational force and perform another iteration. Thus, 
we access the full gravity model at most twice per node while using a number 
of nodes that is substantially lower than in traditional methods. 

Next, let us write the orbit determination problem in a form that conforms 
with the algorithm in Section 14.11 Effectively, we make use of the fact that 
system (134"]) is of the second order. We define the six component vector 



u(t) 



' r(t) ' 




' r(t) ' 


/(*). 




v(f)_ 



where r'(t) = v(£) is the velocity, and the matrix 




where I is 3 x 3 identity matrix. We have 



d 
dt 



' r(t) " 




" r(t) " 







= L 


+ 




v(f) 




v(t) 




-G^ (r(f)) 



(35) 



and the orbit determination problem is now given by (TIB]) with appropriate 
forces as follow from f l3"5"j) . Using f)16p accelerates convergence of the fixed point 
iteration in our scheme. 



4-3 Example 



We present an example of using our method. An extensive study of the method 
for applications in astrodynamics may be found in [3] and here we simply 
demonstrate that our scheme allows large time steps and requires relatively few 
evaluations of the full gravity model. As the cost of evaluating the full (high- 
degree) gravity model is substantial, this results in significant computational 
savings. 

As an example, we simulate an orbit with initial condition 





(r \ 

Xq 


^2284.060 


r\t=o = 


Uo = 


= 6275.400 




w 


v 4.431 
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and 



-5.947 



\ 



dr 

dt \t=o 



2.164 



(km/s), 



V 







and propagate it for 86,000 seconds (approximately 1 day). We use 22 time 
intervals and, on each interval, quadratures with 74 nodes. Hence, on average, 
this corresponds to time distance between nodes of approximately 53 seconds. 
For the full gravitational model we use a 70 degree spherical harmonics model 
WGS84 [5j. 

We compute the orbit trajectory using the algorithm from Sect ion I4TT1 adopted 
to the problem of orbit propagation as described in Section 14.21 We compare 
the result with a trajectory generated by an 8th-order Gauss-Jackson integra- 
tion scheme with one second time step. 

Achieving an error of less than 5 cm at the final time, we need 6512 evaluations 
of the reduced (3-term) gravitational model, and 3256 evaluations of the full 
gravitational model. 



5 Conclusions 



We have constructed an implicit, symplectic integrator that has speed com- 
parable to explicit multistep integrators currently used for orbit computation. 
The key difference with the traditional IRK method is that our scheme uses 
quadratures for band-limited exponentials rather than the traditional Gaus- 
sian quadratures constructed for the orthogonal Legendre polynomials. The 
nodes of quadratures for band-limited exponentials do not concentrate exces- 
sively towards the end points of an interval thus removing a practical limit on 
the number of nodes used within each time interval. 



6 Appendix 

In both approaches described below we use nodes of generalized Gaussian 
quadratures for exponentials {ti} 1=1 constructed in [1] (see Lemma [1]). Some 
of the steps may require extended precision to yield accurate results. 
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6.1 Computing integration matrix using exact PSWFs 



In this approach we assume that the solutions ipj(x) and the eigenvalues \j 
satisfying 

(F c ^) (x) = £ e^(y)dy = X^(x), (36) 

where F c is defined in (j2J), are available. We use Q and the matrix of values 
of PSWFs at the nodes, ipj(ri), to compute coefficients a^j, so that we have 

M-l 

Rt(r)= J2^(r), k = l,...M. 

3=0 

We then compute weights using ([2"2l . 

r i M-l i M-l 

w k = R c k (x)dx = a kj / *Pj(x)dx = a kj^j{0). 
Next we define 

M-i rx M-l 
£ a i3 / ^{s)ds = 

3=0 J 1 j=0 



rx M-l - x M-l 

Ffix) = / R[{s)ds = a U / rf{s)ds = £ aj^^a;), 



where 

* c j (x) = f_r 3 {s)ds. (37) 
In order to compute the integration matrix (|23|) . we need to evaluate 

r i M-l i M-l 

w k S k i= / F l c (x)R c k (x)dx = a ij a kf / ^ c j (x)iljj,(x)dx = "//> /.//,/•• 
1 i,i'=o i,i'=Q 

where 

= J 1 ^ c j (x)rf,(x)dx = j 1 ^ c j (x)—^ c jl (x)dx. (38) 

We have 

Proposition 8 If j and j' are both even, then 

1^ = 1^ = ^X^(0)^(0). (39) 
If j and j' are both odd, then 

I jf = 0. (40) 

If j is even and j' is odd, then 

I33' = -h>v ( 41 ) 
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/ V f 1 ,cf \^f(y)^ 



n icXi j-i 



(42) 



and 



j ' j = ic\~> I ' 1 ^j( y >^~ dy ~ 2 ^K°) 



-i 



2/ 



dy + ic^^Xf / i) c Jy)dy 
o y Jo . 



(43) 



We use (|42]) if \\,\ < |Aj-|. (115)) otherwise. 



PROOF. Integrating ( )38l) by parts, we obtain 

J*, + J yi = ^(1)^,(1) - ^(-l)^(-l) = XjWffltyiP) (44) 
and, since ^ (0) = if j is odd (due to parity of PSWFs), arrive at (1401) and 

m 



Using ( 157)) and ( 155)) . we have 
1 f 1 f f x 



Aj J-i \J-i 



e**-cfa ) ^(y)^ = Y . / : WWv, 

Aj j — i % cy 



and, thus, 



1 f 1 



'"—^J-xiJ-i 



1 e tcyx _ e ~icy 



icy 



rf(y)dy 



ip^,(x)dx 



(45) 



It follows from ( 145)) that if j is even and j' is odd (so that ipy(0) 
obtain (T42]) and 



0), we 



A, / 



y 



' 1 ' ^%)^^%-^ c (0)/ n^-e-^dy 



^%{y) 



i 2/ 



Introducing 



we have 



so that 



u (x) 



j \y ) -icyx 



i y 



dy, 



u\x) = -ic / e "^2/ = —ic\jfif)y(x) 



u(x) = u(a) — icXji I ip^,(s)ds. 

Ja J 

Setting x = 1 and a = 0, we obtain 
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-l y J-i y Jo J 

<hi — 1<-\ > I i , 



2 / -2 dy-icX r / ipUy)dy 

Jo y Jo 



and arrive at (H3 



Computing integration matrix using approximate PSWFs 

If the interpolating basis for band-limited functions is defined via (I10p . then 
the coefficients r k i are obtained using 

M 

5km = Rk(T m ) = Y, r kie iCTlTm (46) 
1=1 

by inverting the matrix E = {e lCTlTm } l m=1 M . We have 

= / R k (s)ds = J2r k i 

J — 1 ; 1 



M_ e icT t x _ e -icri 



l=l 



ICTi 



and compute 



WkSki = J ^Fi(x)R k (x)dx 

r kJ r lf / eiCTjX : dx 

i'=i,...M 1 %CT i' 



where 



3,3 

2J r kj r lj' G jj', 
j,j'=l,...M 



_ sine (c (tj +Tf)) - e iC V s fnc (ctj) 

(jrjji — 2 

ICTji 



Thus, we have 

w k S kl = (E-'GE- 1 
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